home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Pascal Super Library
/
Pascal Super Library (CW International)(1997).bin
/
MATH
/
NS_ROOT
/
NS_ROOTS.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1992-05-11
|
12KB
|
429 lines
{ $A+,B-,D-,F+,I-,L+,N-,R-,S-,V-}
{ $IFOPT N-}
{ $E-}
{ $ELSE}
{ $E+}
{ $ENDIF}
Unit NS_Roots;
{*****************************************************************************}
{
{ TNumMethods
{ +-TRoots
{
{*****************************************************************************}
Interface {===================================================================}
Uses
WinTypes,
WinProcs,
WinDos,
WObjects,
Strings,
StdHdr;
Type
TFloatingPoint = Real; { Format for all floating point numbers }
PNumMethods = ^TNumMethods; {+++++++++++++++++++++++++++++++++++++++++++++++}
TNumMethods = Object(TObject)
ErrorCode : Byte; { Holds most recient error code }
Constructor Init;
Function Min(X1, X2 : TFloatingPoint) : TFloatingPoint; Virtual;
Destructor Done; Virtual;
End; {Object}
PRoots = ^TRoots; {+++++++++++++++++++++++++++++++++++++++++++++++++++++}
TRoots = Object(TNumMethods)
X1 : TFloatingPoint; { Initial lower bracket value }
X2 : TFloatingPoint; { Initial upper bracket value }
Tolerance : TFloatingPoint; { Refine root until + or - Tolerance }
Iter : LongInt; { Iterations performed to converge }
IterMax : LongInt; { Maximum number of iterations allowed (MUST be a positive number!) }
RootValue : TFloatingPoint; { Value at the calculated root }
Constructor Init(vX1, vX2, vTolerance : TFloatingPoint; vIterMax : LongInt);
Function BrentRoots : TFloatingPoint; Virtual;
Function BisectionRoots : TFloatingPoint; Virtual;
Function NewtonRoots : TFloatingPoint; Virtual;
Function fx(X : TFloatingPoint) : TFloatingPoint; Virtual;
Function fxPrime(X : TFloatingPoint) : TFloatingPoint; Virtual;
Destructor Done; Virtual;
Private
OrigX1 : TFloatingPoint;
OrigX2 : TFloatingPoint;
Function Bracketed(vX1, vX2: TFloatingPoint): BOOLEAN;
End; {Object}
Const
nearzero : TFloatingPoint = 1.0e-20;
Implementation {==============================================================}
{++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
{+++ TNumMethods ++++++++++++++++++++++++++++++++++++++++++ TNumMethods +++}
{++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
Constructor TNumMethods.Init;
Begin
TObject.Init;
End;
{EndConstructor}
Function TNumMethods.Min;
{**************************************************************************}
{ Calculates the minimum of two numbers }
{**************************************************************************}
BEGIN
IF (X1 < X2) THEN
Min := X1
ELSE
Min := X2;
{EndIf}
END;
{EndFunction}
Destructor TNumMethods.Done;
Begin
End;
{End Destructor}
{+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
{+++ TRoots +++++++++++++++++++++++++++++++++++++++++++++++++++ TRoots +++}
{+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++}
Constructor TRoots.Init;
{**************************************************************************}
Begin
{ NumMethods.Init;}
X1 := vX1;
OrigX1 := vX1;
X2 := vX2;
OrigX2 := vX2;
Tolerance := vTolerance;
{--------------------------------------------------------------------}
{ We don't want negative itterations! If vIterMax is less than 1 }
{ set IterMax to the largest value a LongInt will hold. }
{--------------------------------------------------------------------}
If (vIterMax < 1) Then
IterMax := 2147483647
Else
IterMax := vIterMax;
{EndIf}
End;
{EndConstructor}
Function TRoots.fx(X: TFloatingPoint): TFloatingPoint;
{**************************************************************************}
BEGIN
Abstract;
END;
{EndFunction}
Function TRoots.fxPrime(X: TFloatingPoint): TFloatingPoint;
{**************************************************************************}
BEGIN
Abstract;
END;
{EndFunction}
Function TRoots.BrentRoots;
{**************************************************************************}
CONST
Finished : Boolean = False;
FPP = 1.0e-11; { Machine floating point precision }
VAR
result, CC, DD, EE, fa, fb, fc, Tol1, PP, QQ, RR, SS, XM: TFloatingPoint;
BEGIN
Iter := 0;
X1 := OrigX1;
X2 := OrigX2;
ErrorCode := 0;
fa := fx(X1);
fb := fx(X2);
IF (NOT(Bracketed(fa,fb))) THEN
ErrorCode := 1
ELSE
BEGIN
fc := fb;
REPEAT
IF (NOT(Bracketed(fc,fb))) THEN
Begin
CC := X1;
fc := fa;
DD := X2 - X1;
EE := DD;
End;
{EndIf}
IF (ABS(fc) < ABS(fb)) THEN
BEGIN
{---------------------------------------------------}
{ Rename X1, X2, CC and adjust bounding interval DD }
{---------------------------------------------------}
X1 := X2;
X2 := CC;
CC := X1;
fa := fb;
fb := fc;
fc := fa;
END;
{EndIf}
Tol1 := 2.0 * FPP * ABS(X2) + 0.5 * Tolerance; { Convergence check }
XM := 0.5 * (CC-X2);
IF ((ABS(XM) <= Tol1) OR (ABS(fa) < nearzero)) THEN
BEGIN
result := X2;
Finished := TRUE;
RootValue := fx(result);
END
ELSE
BEGIN
IF ((ABS(EE) >= Tol1) AND (ABS(fa) > ABS(fb))) THEN
BEGIN
SS := fb/ fa;
IF (ABS(X1 - CC) < nearzero) THEN
BEGIN
PP := 2.0 * XM * SS;
QQ := 1.0 - SS;
END
ELSE
BEGIN
QQ := fa/fc;
RR := fb /fc;
PP := SS * (2.0 * XM * QQ * (QQ - RR) - (X2-X1) * (RR - 1.0));
QQ := (QQ - 1.0) * (RR - 1.0) * (SS - 1.0);
END;
{EndIf}
IF (PP > nearzero) THEN
QQ := -QQ;
{EndIf}
PP := ABS(PP);
IF ((2.0 * PP) < Min(3.0*XM *QQ-ABS(Tol1 * QQ), ABS(EE * QQ))) THEN
Begin
EE := DD;
DD := PP/QQ;
End
ELSE
Begin
DD := XM;
EE := DD;
End
{EndIf}
END
ELSE
BEGIN
DD := XM;
EE := DD;
END;
{EndIf}
X1 := X2;
fa := fb;
IF (ABS(DD) > Tol1) THEN
X2 := X2 + DD
ELSE
IF (xm > 0.0) THEN
X2 := X2 + ABS(Tol1)
ELSE
X2 := X2 - ABS(Tol1);
{EndIf}
{EndIf}
fb := fx(X2);
Inc(Iter);
IF (Iter >= IterMax) THEN
Begin
ErrorCode := 2;
Finished := True;
End
{EndIf}
END;
{EndIf}
UNTIL (Finished);
END;
{EndIf}
BrentRoots := result;
END;
{EndMethod}
Function TRoots.BisectionRoots;
{**************************************************************************}
CONST
FPP = 1.0e-11;
nearzero = 1.0e-20;
Finished : Boolean = False;
VAR
result : TFloatingPoint;
fa, fb, fc : TFloatingPoint;
XC : TFloatingPoint;
BEGIN
Iter := 0;
X1 := OrigX1;
X2 := OrigX2;
ErrorCode := 0;
fa := fx(X1);
fb := fx(X2);
IF (NOT(Bracketed(fa, fb))) THEN
ErrorCode := 1
ELSE
BEGIN
REPEAT
XC := (X1 + X2) / 2.0;
fc := fx(XC);
IF (Bracketed(fa, fc)) THEN
X2 := XC
ELSE
Begin
X1 := XC;
fa := fc;
End;
{EndIf}
IF (ABS(X2 - X1) <= Tolerance) THEN
BEGIN
Finished := TRUE;
result := XC;
RootValue := fx(result);
END
Else
Begin
Inc(Iter);
IF (Iter >= IterMax) THEN
Begin
ErrorCode := 2;
Finished := True;
End;
{EndIf}
End;
{EndIf}
UNTIL (Finished);
END;
{EndIf}
BisectionRoots := result;
END;
{EndMethod}
Function TRoots.NewtonRoots;
{**************************************************************************}
CONST
Finished : BOOLEAN = False;
FPP = 1.0e-11;
nearzero = 1.0e-20;
VAR
result, lastX, fxCalc, FPX, XP : TFloatingPoint;
BEGIN
X1 := OrigX1;
X2 := OrigX2;
ErrorCode := 0;
fxCalc := fx(X1);
Iter := 0;
REPEAT
IF (abs(X1) < nearzero) THEN
Begin
ErrorCode := 1;
Finished := True;
End
ELSE
BEGIN
FPX := fxPrime(X1);
IF (abs(FPX) < nearzero) THEN
Begin
ErrorCode := 2;
Finished := True;
End
ELSE
BEGIN
lastX := X1;
X1 := X1 - fxCalc/FPX;
fxCalc := fx(X1);
IF (ABS(X1 - lastX) <= Tolerance) THEN
BEGIN
Finished := TRUE;
result := X1;
RootValue := fx(X1);
END
Else
Begin
Inc(Iter);
IF (Iter >= IterMax) THEN
Begin
ErrorCode := 3;
Finished := True;
End
{EndIf}
End
{EndIf}
END
{EndIf}
END;
{EndIf}
UNTIL (Finished);
NewtonRoots := result;
END;
{EndMethod}
Function TRoots.Bracketed;
{**************************************************************************}
{ Return True if vX1 * vX2 is negative. }
{**************************************************************************}
BEGIN
IF ((vX1 > 0.0) AND (vX2 > 0.0)) OR ((vX1 < 0.0) AND (vX2 < 0.0)) THEN
Bracketed := FALSE
ELSE
Bracketed := TRUE;
{EndIf}
END;
{EndFunction}
Destructor TRoots.Done;
{**************************************************************************}
Begin
TNumMethods.Done;
End;
{EndDestructor}
End.